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_C3 Abstract 

^ This paper addresses the issue of structure-preserving discretization of open distri- 

buted-parameter systems with Hamiltonian dynamics. Employing the formahsm of 
discrete exterior calculus, we introduce a simplicial Dirac structure as a discrete 
analogue of the Stokes-Dirac structure and demonstrate that it provides a natural 
^ framework for deriving finite-dimensional port-Hamiltonian systems that emulate their 

CO infinite-dimensional counterparts. The spatial domain, in the continuous theory rep- 

resented by a finite-dimensional smooth manifold with boundary, is replaced by a 
homological manifold-like simplicial complex and its augmented circumcentric dual. 
The smooth differential forms, in discrete setting, are mirrored by cochains on the 
primal and dual complexes, while the discrete exterior derivative is defined to be the 
^ coboundary operator. This approach of discrete differential geometry, rather than 

I discretizing the partial differential equations, allows to first discretizc the underlying 

Stokes-Dirac structure and then to impose the corresponding finite-dimensional port- 
. ^ Hamiltonian dynamics. In this manner, a number of important intrinsically topological 

and geometrical properties of the system are preserved. 



1 Introduction 



The purpose of this paper is to propose a sound geometric framework for structure- 
preserving discretization of distributed-parameter port-Hamiltonian systems. Our ap- 
proach to time-continuous spatially-discrete port-Hamiltonian theory is based on discrete 
exterior geometry and as such proceeds ab initio by mirroring the continuous setting. The 
theory is not merely tied to the goal of discretization but rather aims to offer a sound 
and consistent framework for defining port-Hamiltonian dynamics on a discrete manifold 
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which is usually, but not necessarily, obtained by discretization of a smooth Riemannian 
manifold. 

The underlying structure of open distributed-parameter dynamical systems considered 
in this paper is a Stokes-Dirac structure [26] and as such is being defined on a certain 
space of differential forms on a smooth finite-dimensional orientable, usually Riemannian, 
manifold with a boundary. The Stokcs-Dirac structure generalizes the framework of the 
Poisson and symplectic structures by providing a theoretical account that permits the 
inclusion of varying boundary variables in the boundary problem for partial differential 
equations. Prom an interconnection and control viewpoint, such a treatment of boundary 
conditions is essential for the incorporation of energy exchange through the boundary, 
since in many applications the interconnection with the environment takes place precisely 
through the boundary. The same arguments apply to the finite-dimensional approxima- 
tions of complex distributed-parameter systems. Por numerical integration, simulation 
and control synthesis, it is of paramount interest to have finite approximations that can 
be interconnected to one another or via the boundary coupled to other systems, be they 
finite- or infinite-dimensional. 

Most of the numerical algorithms for spatial discretization of distributed-parameter 
systems, primarily finite difference and finite element methods, fail to capture the intrinsic 
system structures and properties, such as symplecticity, conservation of momenta and 
energy, as well as differential gauge symmetry. Purthermore, some important results, 
including the Stokes theorem, fail to apply numerically and thus lead to spurious results. 
This loss of fidelity to preserve some inherent topological and geometric structures of the 
continuous models motivates a more geometry based approach. 

The discrete approach to geometry goes back to Whitney, who in [38] introduced an 
isomorphism between simplicial and de Rham cohomology. More recent antecedents can 
be found, for instance, in [30], and also in the computational electromagnctism literature 
[3, 4, 14]. Por a comprehensive historical summary we refer to the thesis [15] and refer- 
ences therein. The literature, however, seems mostly focused on discretization of systems 
with infinite spatial domains, boundaryless manifolds, and systems with zero boundary 
conditions. In this paper, we augment the definition of the dual cell complex in order to 
allow nonzero energy fiow through boundary. 

A notable previous attempt to resolve the problem of structure-preserving discretiza- 
tion of port-Hamiltonian systems is [12], where the authors employ the mixed finite element 
method. Their treatment is restricted to the one-dimensional telegraph equation and the 
two-dimensional wave equation. Although it is hinted that the same methodology applies 
in higher dimensions and to the other distributed-parameter systems, the results are not 
clear. It is worth noting that the choice of the basis functions can have dramatic conse- 
quences on the numerical performance of the mixed finite element method; as the mesh 
is being refined, it easily may lead to an ill-conditioned finite-dimensional linear system 
[2] . The other undertaking on discretization of port-Hamiltonian systems can be found in 
[27, 28] , but the treatment is purely topological and is more akin to the graph-theoretical 
formulation of conservation laws. Purthermore, the authors in [27, 28] do not introduce 
a discrete analogue of the Stokes theorem and the entire approach is tied to the goal of 
preserving passivity. 

Our approach is that of discrete exterior calculus [8, 9, 15, 33], which has previously 
been applied to variational problems naturally arising in mechanics and electromagnetism. 
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These problems stem from a Lagrangian, rather than Hamiltonian, modehng perspec- 
tive and as such they conform to a muhisymplectic structure [13, 22, 23, 35], rather 
than the Stokes-Dirac structure. A crucial ingredient for the numeric integration is the 
asynchronous variational integrator for spatio-temporally discretized problems, whereas 
our approach spatially discretizes the Stokes-Dirac structure and allows imposing time- 
continuous spatially discrete dynamics. This apparent discrepancy between multisym- 
plectic and the Stokes-Dirac structure-preserving discretization could be elevated by, for 
instance, defining Stokes-Dirac structure on a pseudo-Riemannian manifold to insure a 
treatment of space and time on equal footing, whilst keeping nonzero exchange through 
the boundary. 

Contribution and outline of the paper. We begin by recalling the definition of the 
Stokes-Dirac structure and port-Hamiltonian systems. In order to make this paper as self- 
contained as possible for a variety of readers, we present a brief overview of the elementary 
discrete exterior geometry needed to define a discretized Stokes-Dirac structure and im- 
pose appropriate port-Hamiltonian dynamics. The third section is a brief summary of the 
essential definitions and results in discrete exterior calculus as developed in [8, 9, 15]. The 
contribution of this paper in this regard is a proper treatment of the boundary of the dual 
cell complex. Namely, in order to allow the inclusion of nonzero boundary conditions on 
the dual cell complex, we off'er a definition of the dual boundary operator that differs from 
the standard one. Such a construction leads to a discrete analogue of the integration by 
parts formula, which is a crucial ingredient in establishing a discrete Stokes-Dirac structure 
on a primal simplicial complex and its circumcentric dual. The main result is presented in 
Section 4, where we introduce the notion of simplicial Dirac structures on a primal-dual 
cell complex, and in the following section define port-Hamiltonian systems with respect 
to these structures. In Section 6 we give a matrix representation for the simplicial Dirac 
structures and linear port-Hamiltonian systems, for which we also establish bounds for the 
energy of discretization errors. Finally, wc demonstrate how the simplicial Dirac struc- 
tures relate to some spatially discretized distributed-parameter systems with boundary 
variables: Maxwell's equations on a bounded domain, a two-dimensional wave equation, 
and the telegraph equations. While the focus of this paper is not implementation of dis- 
crete exterior calculus in discretization of port-Hamiltonian systems, wc have, nonetheless, 
taken some preliminary numerical investigations. We demonstrate the application of the 
developed machinery using the example of the telegraph equations. 
Some preliminary results of this paper have been reported in [32] . 

2 Dirac Structures and Port-Hamiltonian Dynamics 

Dirac structures were originally developed in [5, 6, 11] as a generalization of symplectic 
and Poisson structures. The formalism of Dirac structure was employed as the geometric 
notion underpinning generalized power-conserving interconnections and thus allowing the 
Hamiltonian formulation of interconnected and constrained dynamical systems. 

A constant Dirac structure can be defined as follows. Let J^, £, and L be linear spaces. 
Given a / G J" and an e G the pairing will be denoted by (e|/) G L. By symmetrizing 
the pairing, we obtain a symmetric bilinear form ((, )) : F x £ ^ L defined by 

(((/i,ei),(/2,e2))) = (ei|/2) + (e2|/i). 
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Definition 2.1. A Dirac structure is a linear subspace T> d F x £ such that V = , 
with _L standing for the orthogonal complement with respect to the bilinear form ((, )) . 

It immediately follows that for any (/, e) G P 

0=(((/,e),(/,e))) = 2(e|/). 

Interpreting (/, e) as a pair of power variables, the condition (/, e) G I? implies power- 
conservation (e|/) = 0, and as such is terminus a quo for the geometric formulation of 
port-Hamiltonian systems. 

Much is known about finite-dimensional Dirac structures and their role in physics; 
however, hitherto there is no complete theory of Dirac structures for field theories. An 
initial contribution in this direction is made in the paper [26] , where the authors introduce 
a notion of the Stokes-Dirac structure. This infinite-dimensional Dirac structure lays 
down the foundation for port-Hamiltonian formulation of a class of distributcd-parameter 
systems with boundary energy flow. In this section, we provide a very brief overview of 
the Stokes-Dirac structure [26]. 

Throughout this paper, let M be an oriented n-dimensional smooth manifold with a 
smooth (n — l)-dimcnsional boundary dM endowed with the induced orientation, repre- 
senting the space of spatial variables. By ^^^{M), k = 0,1,..., n, denote the space of 
exterior A;-forms on M, and by Qf^{dM), = 0, 1, . . . , n — 1, the space of fc-forms on dM. 
A natural non-degenerative pairing between a G Q,^{M) and (3 G Q'^'~^{M) is given by 
(/3|a) = (3 A a. Likewise, the pairing on the boundary dM between a G r2^(5M) and 
/3 G f^^-^-^SM) is given by (/?|a) = Jq^ P A a. 

For any pair p, q of positive integers satisfying p + q = n + 1, define the flow and effort 
linear spaces by 

Sp,q=n''-P{M) X J)"-«(M) X n"-«(5M) . 
The bilinear form on the product space J-'p^q x £p^q is given by 

iiifp^fq ' fhi^p^ ^b)i ifp^fq > fb^^p^ ^&))) = 
V ' V ' 

&^p,q ^£p,q (2.1) 

/ [elAfl + e\Afl + elAfl+elAfl)+ I [el A + el A fl) . 

JM JdM 

Theorem 2.1. Given linear spaces Fp^q and £p^q, and bilinear form ((,)), define the fol- 
lowing linear subspace P of Fp^q x £p^q 

P ~ { fq^fbjGp,eq, efe) G -^p.q ^ ^p^q \ 



fq J \d J \eq J ' (2.2) 

Cb J \0 -(-!)"-« J \ eq\9M 



where d is the exterior derivative and \gM stands for a trace on the boundary dM. Then 
T> = T>^ , that is, V is a Dirac structure. 
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Remark 2.1. Although the differential operator in (2.2), in the presence of nonzero bound- 
ary conditions, is not skew- symmetric, it is possible to associate a pseudo-Poisson structure 
to the Stokes-Dirac structure [26]. In the absence of algebraic constraints, the Stokes-Dirac 
structure specializes to a Poisson structure [1 1 ], and as such it can be derived through sym- 
metry reduction from a canonical Dirac structure on the phase space [36]. Whether this 
reduction can be done for the Stokes-Dirac structure on a manifold with boundary remains 
an important open problem. 

In order to define Hamiltonian dynamics, consider a Hamiltonian density Ti : J7^(M) x 
Q'?(M) — > 0"(M) resulting with the Hamiltonian H = fj^H G R. Now, consider a 
time function t i-)- {ap{t) , aq{t)) G Qp{M) x J^«(M), t G M, and the Hamiltonian t 
H {ap{t) , aq{t)) evaluated along this trajectory, then at any t 



where {SpH,6qH) e Jl" ^(M) x *(M) are the (partial) variational derivatives of H at 



Setting the flows fp = — g^, fg = —-Qf and the efforts Cp = 5pH, Cq = 5qH, the 
distributed-parameter port-Hamiltonian system is defined by the relation 



For such a system, it straightaway follows that ^ = Jqj^ Cb A ft, expressing the fact that 
the system is lossless. In other words, the increase in the energy of the system is equal to 
the power supplied to the system through the boundary dM. 

3 Fundaments of Discrete Exterior Calculus 

The discrete manifolds we employ are oriented manifold-like simplicial complexes and their 
circumccntric duals. Typically, these manifolds are simplicial approximations of smooth 
manifolds. Familiar examples are meshes of triangles embedded in and tetrahedra 
obtained by tetrahedrization of a 3-dimensional manifold. There are many ways to obtain 
such complexes; however, we do not address the issue of discretization and embedding. 

As said, we proceed ab initio and mostly our treatment is purely formal, that is without 
proofs that the discrete objects converge to the continuous ones, though we briefly address 
the issue of convergence in Section 6.3. By construction of discrete exterior calculus, a 
number of important geometric structures are preserved and propositions like the Stokes 
theorem are true by definition. The basic building blocks of discrete exterior geometry are 
discrete chains and cochains, and their geometric duals. The former are simplices and the 
latter are discrete differential forms related one to another by bilinear pairing that can be 
understood as the evaluation of a cochain on an appropriate simplex, and as such parallels 
the integration in the continuous setting. 

In discrete exterior calculus, a dual mesh is instrumental for defining the diagonal 
Hodge star. In the paper at hand, the geometric duality is a crucial ingredient in estab- 
lishing a bijective relationship between the flow and effort spaces, as well as for construction 
of nondegenerate discrete analogues of the bilinear form (2.1). 
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This section, with some modification concerning the treatment of the boundary of the 
dual cell complex, is a brief summary of the essential definitions and results in discrete 
exterior calculus as developed in [8, 9, 15]. As therein, first we define discrete differential 
forms, the discrete exterior derivative, the codifferential operator, the Hodge star, and the 
discrete wedge product. For more information on the construction of the other discrete 
objects such as vector fields, a discrete Lie derivative, and discrete musical operators, we 
refer the reader to [15]. Construction of all these discrete objects is in a way simpler than 
their continuous counterparts since we require only a local metric, ergo the machinery of 
the Riemannian geometry is not demanded. With the exception of the treatment of the 
notions related the boundary of the dual cell complex, a good part of this section is a 
recollection of the basic concepts and results of algebraic topology [16, 25]. 

3.1 Simplicial complexes and their circumcentric duals 

Definition 3.1. A k-simplex is the convex span ofk + 1 geometrically independent points, 



The points vo, ■ ■ ■ ,Vk are called the vertices of the simplex, and the number k is called the 

dimension of the simplex. Any simplex spanned by a (proper) subset of {vq, . . . ,Vk} is 
called a (proper) face of . If is a proper face of a^, we denote this by ~< . 

As an illustration, consider four non-collinear points vq, V2, and in R"'. Each of 
these points individually is 0-simplex with an orientation dictated by the choice of a sign. 
An example of a 1-simplex is a line segment [t'o,t'i] oriented from vq to vi. The triangle 
[vq, ^1,^2] is an example of 2-simplex oriented in counterclockwise direction. Similarly, the 
tetrahedron [vq,vi,V2,vz\ is a 3-simplex. 

Definition 3.2. A simplicial complex K in is a collection of simplices in M.^ , such 
that: 

(1) Every face of a simplex of K is in K. 

(2) The intersection of any two simplices of K is a, face of each of them. 

The dimension n of the highest dimension simplex in K is the dimension of K. 

The above given definition of a simplicial complex is more general than needed for the 
purposes of exterior calculus. Since the discrete theory employed in this paper mirrors the 
continuous framework, we restrict our considerations to manifold-like simplicial complexes 



Definition 3.3. A simplicial complex K of dimension n is a manifold-like simplicial 
complex if the underlying space is a polytope \K\. In such a complex all simplices of 
dimension A; = 0, . . . , n — 1 must be a face of some simplex of dimension n in the complex. 

Introducing these simplicial meshes has an added advantage of allowing a simple and 
intuitive definition of orientability of simplicial complexes [15]. 




[15]. 
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Definition 3.4. An n- dimensional simplicial complex K is an oriented manifold-like sim- 
plicial complex if the n-simplices that share a common (n — l)-face have the same orien- 
tation and all the simplices of lower dimensions are individually oriented. 

Henceforth in this paper, we shall work with manifold-like simplicial complexes. When 
no confusion can arise, we address these objects simply as simplicial complexes. 

An essential constituent of discrete exterior calculus is the dual complex of a manifold- 
like simplicial complex. The most popular notions of duality are barycentric and cir- 
cumcentric, also known as Voronoi, duality. Following the standard approach of discrete 
exterior calculus, in this paper we employ the latter. 

The circumcenter of a fc-simplex is given by the center of the A;-circumsphere, where 
the fc-circumsphere is the unique fc-sphere that has all A; + 1 vertices of on its surface. 
That is, the circumcenter is the unique point in the fc-dimensional affine space that contains 
the fc-simplex that is equidistant from all the fc + 1 nodes of the simplex. We denote the 
circumcenter of a simplex by c{a^). 

If the circumcenter of a simplex lies in its interior we call it a well-centered simplex. 
For instance, a triangle with all acute angles is a well-centered 2-simplex. A simplicial 
complex K whose all simplices of all dimensions are well-centered is called a well-centered 
simplicial complex and its dual obtained by circumcentric subdivision is also a simplicial 
complex denoted by csdK and its elements by (T°,...,(t". Throughout this paper, we 
adopt the convention that all symbols related to the dual (simplicial and cell) complex are 
labeled by a caret. The underlying spaces \K\ and |csdi^| are the same. The simplicial 
complex csd K consists of all simplices of the form [c((Ti), . . . , c((7fc)] for = 1, . . . , n, where 
o"! ^ (72 ^ ... ^ CTfc, meaning ui is a proper face of Uj for all i < j. 

A circumcentric dual cell complex (block complex in terminology of [25]) is obtained 
by aggregation of certain simplices of csd K. Let K he a well-centered simplicial complex 
of dimension n and let a'' be one of its simplices. By D{aP) we denote the union of all 
open simplices of csdK of which c{a'^) is the final vertex; this cell is the dual cell to a'^. 
The closure of the dual cell of a'^ is D{a^). The collection of all dual cells is a cell complex 
denoted by D{K) with closure D{K). 

To illustrate the duality, consider the 2-dimensional simplicial complex pictured in 
Figure 1. The dual cell of the vertex Vr is the topological interior of the Voronoi region 
around it as shown shaded in the figure. This dual cell is comprised of the vertex w^, 
the interior of the open edges emanating from Vr, and interiors of the all dual simplices 
containing Vr- The dual cell of any 2-simplex consists of its circumcenter alone. The dual 
cell of an edge consists of its circumcenter and two open edges joining this circumcenter 
to the circumccntcrs of two triangles having the primal edge as a face. The dual of a 
boundary edge has only one half-edge since there is only one triangle adjacent to that 
boundary edge. Note that if the complex is not flat, then the dual edge will not be a 
straight line. 

Remark 3.1. A triangulation of a compact n-dimensional Riemannian manifold M re- 
sults in an n-dimensional simplicial complex K. Intuitively, the simplices are glued to the 

manifold M in such a way that they form a 'curved' manifold-like simplicial complex. It is 
worth noticing that in practical applications, the smooth manifold sometimes is unknown 
and can only be sampled by physical measurements. In such situations, it makes sense to 
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Figure 1: A 2-dimensional simplicial complex K subdivided into the circumcentric simpli- 
cial complex csdK indicated by dotted lines. The dual cells displayed are shaded. 

model the spatial domain as inherently discrete. This is where discrete port-Hamiltonian 
theory in the framework of discrete exterior calculus stands in its own right. 

3.2 Chains and cochains 

The discrete analogue of a smooth k-foim is a fc-cochain, a certain type of a function, on 
a /c-chain representing a formal sum of simplices. The role of integration in the discrete 
theory is replaced by (simple) evaluation of a discrete form on a chain. The discrete exterior 
derivative is defined by duality to the boundary operator, rendering the Stokes theorem 
true by definition. Parallel to the smooth case, the discrete exterior wedge product pairs 
lower degree forms into a higher degree one. 

Definition 3.5. Let K be a simplicial complex. We denote the free Abelian group gener- 
ated by a basis consisting of oriented k-simplices by CkiK^Tj). This is the space of finite 
formal sums of the k-simplices with coefficients in Z. Elements of Ck{K;Z) are called 
k-chains. 

Definition 3.6. A primal discrete k-form a is a homomorphism from the chain group 
Ck{K; Z) to the additive group M. Thus, a discrete k-form is an element ofHom{Ck{K),M.), 
the space of cochains. This space becomes an Abelian group if we add two homomorphisms 
by adding their values in M. The standard notation for Hom(Cfc(ir), M) in algebraic topol- 
ogy is C^{K;M.); however, like in [8, 9, 15] we shall also employ the notation ^1^{K) for 
this space as a reminder that this is the space of discrete k- forms on the simplicial complex 
K. Thus, 

^d{K) := C^{K;M.) = Hom(Cfc(K),M) . 
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Given a /c-chain aicl, ai G Z, and a discrete A;-form a, we have 




and for two discrete A;-forms a, /3 G ^^{K) and a A;-chain c G Ck{K; Z), 

(a + /3)(c) = a(c) + /?(c). 

The natural pairing of a A;-form a and a fe-chain c is defined as the bilinear pairing 
{a, c) = a(c). 

As previously pointed out, a differential A;-form can be thought of as a linear 
functional that assigns a real number to each oriented cell G K. In order to understand 
the process of discretization of the continuous problem consider a smooth /c-form / G 
Jl'^diirl). The discrete counterpart of / on a fc-simplex ^ K \s & discrete form 
defined as a{a^) := J^k /• 

Definition 3.7. The boundary operator dk ■ CkiK]Z) — )• Ck-i{K;Z) is a homomorphism 
defined by its action on a simplex a'^ = [vq, . . . ,Vk], 

k 

dkcr'' = dk{[vo,...,Vk]) = ^{-iy[vo, . . . ,Vi, . . . ,Vk] , 

where [vq, . . . ,Vi, . . . ,Vk] is the {k — l)-simplex obtained by omitting the vertex vi. Note 
that dk o dk+i = 0. 

Definition 3.8. On a simplicial complex of dimension n, a chain complex is a collection 
of chain groups and homomorphisms dk, such that 

0—^CniK) ^ Cn-l ^ ■ ■ ■ ^ Ck{K) Co{K) ^ , 

and dk o dk+i = 0. 

Definition 3.9. The discrete exterior derivative d : — t- Q^'^^(K) is defined by 
duality to the boundary operator dk+i '■ Ck+i{K;Z) — t- Ck{K;Z), with respect to the natural 
pairing between discrete forms and chains. For a discrete form G Q,^(K) and a chain 
Ck+i G Cfe+i {K; Z) we define d by 

(da'',CA;+i) = (a'',afc+iCfc+i) . 

The discrete exterior derivative is the coboundary operator from algebraic topology [25] 
and as such it induces the cochain complex 

o^n-,{K) ^ ny' ^ • • . ^ j^O(K) ^0, 

where d o d = 0. 
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Such as in the continuous theory, we drop the index of the boundary operator when 
its dimension is clear from the context. The discrete exterior derivative d is constructed 
in such a manner that the Stokes theorem is satisfied by definition. This means, given a 
{k + l)-chain c and a discrete A;-form a, the discrete Stokes theorem states that 

(da, c) = {a, dc) . 

Consider a /c-chain ^^aici, G Z, Cj G Cfc(i^;Z), and [k — l)-form a G Q^^{K;'L). 
By hnearity of the chain-cochain pairing, the discrete Stokes theorem can be stated as 



da, ^ aiCi ) = ( ^ XI ) " ( XI ) " X ^^"^ 




As in the continuous setting, the discrete wedge product pairs two discrete differential 
forms by building a higher degree form. The primal-primal wedge product inherits some 
important properties of the cup product such as the bilinearity, anticommutativity and 
naturality under pull-back [8, 15]; however, it is in general non-associative and degenerate, 
and thus unsuitable for construction of canonical pairing between the flow and effort space. 
For a definition of nondegcncratc pairing between the flow and effort discrete forms we 
shall use a primal-dual wedge product as will be defined in the subsequent section. 

3.3 Metric-dependent peirt of discrete exterior calculus 

A cellular chain group associated with the dual cell complex D{K), in [25] denoted by 
Dp{K), is the group of formal sums of cells with integer coefficients. Since in D{K) 
the information of dual simpliccs is lost, to retain the bookkeeping information Hirani in 
[15] introduces a duality operator which takes values in the domain group Cp(csdK;Z). 
As will be clear from the subsequent section, this bookkeeping is not indispensable for 
the formulation of the Dirac structure on a simplicial complex; nevertheless, since the 
information of dual simplices might be needed in defining dynamics, we also employ this 
construction. 

In order to explicitly construct the duality on the boundary, in the next definition 
we introduce the boundary star operator. Shortly afterward we shall explain the rational 
behind this construction. 

Definition 3.10. Let K he a well-centered simplicial complex of dimension n. The interior 
circumcentric duality operator -ki : Ck^K;"!,) — > C„_fe(csdi^;Z) 

= X [c{a%c{a^+^l . . . , c{a-) 

and the boundary star operator ~k\) : Ck{dK;7.) — > C„_i_fc(5(csdK);Z) 

where the s^k and s^k ^n-i coefficients ensure that the orientation of the cell [c{a'^), 
c(cr'^+^), . . . , c((t")] and [c{a''),c{a^^^), . . . ,c{a"'~^)] is consistent with the orientation of 
the primal simplex, and the ambient volume forms on K and dK, respectively. 
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Figure 2: The circumcentric dual cell complex -kK of the simplicial complex K given in 
Figure 1. The boundary of -kK is the dual of the boundary of K. Some support volumes 
are shaded. For instance, the support volume of the primal vertex Vr is the area of its 
Voronoi region; also Vj^^^^^j = V^q^^ ,,^,]). 

The subset of chains Cp(csd K] Z) that are equal to the cells of D{K) x D(dK) forms a 
subgroup of Cp(csdK;Z). We denote this subgroup of Cp(csdi^;Z) by Cp(*i^;Z), where 
-kK is its basis set. A cell complex -kK in is a collection of cells in such that: (1) 
there is a partial ordering of cells in -kK, &^ ~< cr', which is read as is a face of (2) 
the intersection of any two cells in krK, is either a face of each of them, or it is empty; (3) 
the boundary of a cell is expressed as a sum of its proper faces. 

Given a simplicial well-centered complex K, we define its interior dual cell complex 
•kiK (block complex in terminology of algebraic topology [25]) as a circumcentric dual 
restricted to \K\. An important property of the the Voronoi duality is that primal and 
dual cells are orthogonal to each other. The boundary dual cell complex *b-f^ is a dual 
to dK. The dual cell complex -kK is defined as -kK = -k\K x -k\^K. A dual mesh ★iiC 
is a dual to K in sense of a graph dual, and the dual of the boundary is equal to the 
boundary of the dual, that is d{-kK) = •k(dK) = -k^^K. This construction of the dual is 
compatible with [33, 18] and as such is very similar to the use of the ghost cells in finite 
volume methods in order to account for the duality relation between the Dirichlet and the 
Neumann boundary conditions. Because of duality, there is a one-to-one correspondence 
between /c-simplices of K and interior (n — /c)-cells of -kK. Likewise, to every A:-simplex on 
dK there is a uniquely associated (n — 1 — A;)-cell on d{-kK). 

In what follows, we shall abuse notation and use the same symbol * for both the 
interior circumcentric and the boundary star operator. The difference, if not clear from 
the exposition, will be delineated by indicating that -ka^ £ d{-kK) when -ka^ is a dual cell 
on the boundary of the dual cell complex -kK. As sets, the set of D{a^) and -ka^ are equal, 
the only difference being in the bookkeeping, since in -ka^ one retains the information 
about the simplices it is built of. Here we do not address the problem of the orientation 
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of dual -kK, for which we direct the reader to [15]. The circumcentric dual cell complex of 
the 2-dimensional simpHcial complex from Figure 1 is pictured in Figure 2. 

An important concept in defining a wedge product between primal and dual cochains 
is the notion of a support volume associated with a given simplex or cell. 

Definition 3.11. The support volume of a simplex is an n-volume given by the convex 
hull of the geometric union of the simplex and its circumcentric dual. This is given by 

V^u = CH(a^★i(7'=)n|i^|, 

where CH((t'', *iCJ*^) is the n-dimensional convex hull generated by cr*^ U^icr*^. The inter- 
section with \K\ is necessary to ensure that the support volume does not extend beyond 
the polytope \K\ which would otherwise occur if K is nonconvex. The support volume of a 
dual cell -kia^ is 

V^.^k = CH(*i(7'=,*i a^) n \K\ = V^k . 

Everything that has been said about the primal chains and cochains can be extended 
to dual cells and dual cochains. We do not elaborate on this since it can be found in 

the literature [15, 9], however, in order to properly account for the behaviours on the 
boundary, we need to adapt the definition of the boundary dual operator as presented in 
[15, 9]. We propose the following definition. 

Definition 3.12. The dual boundary operator dk ■ Ck{-*:iK;Z) -> Cfc_i(*i^;Z) is a homo- 
morphism defined by its action on a dual cell = ★icr""'^ = -ki[vo, . . . , v^-k] , 

da'' = d-ki [vq, . . . , Vn-k] = di -ki [vo, Vn-k] + db*i [vo, . . . , Vn-k] , 

where 

di-ki [vo,...,Vn-k] 
d^-ki [vo,...,Vn-k] 

Note that the dual boundary operator as defined in [15] is equal to di. Hence, the 
dual boundary is not the geometric boundary of a cell, because near the boundary of a 
manifold that would be wrong. As an example consider the complex in Figure 3. The 
dual of the vertex vi is the Voronoi region shown shaded. Its geometric boundary has five 
sides (two half primal edges and three dual edges), whereas the dual boundary according 
to the definition given in [15] consists of just dual edges, i.e. [c([fo, t'l]), c([wo, fi, ^2])], 
[c{[vo,vi,V2]),c{[vi,V3,V2])] and [c([vi, ^3, W2]), c([ui, U3])], ah up to a sign depending on 
the chosen orientation. However, according to the above given definition, the boundary is 
comprised of four edges, three already given plus the boundary edge [c{[vo, vi]), c{[vi,V3])] 
obtained by aggregation of the two dual simplices [c([uo, wi]), f 1] and [fi, c([t;i, t^a])] . This 
construction of the dual boundary ensures a natural pairing between a primal 0-form 
defined on vi and a dual 1-form on [c([fo, Ui]), c([fi, W3])]. The offered definition of the 
dual boundary operator, as will be demonstrated later, is crucial for the inclusion of the 
boundary variables in the discrete setting. 



fj-n—k+lyfj-n—k 
= *b (s^n-kd''-'') . 
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Figure 3: A 2-dimensional simplicial complex taken from Figure 3.3, Section 3.6, [15]. The 
shaded region is a Voronoi dual of the primal vertex vi. The dual boundary, according to 
[15], is not the geometric boundary near the boundary of the manifold. However, in line 
with our construction, the boundary of the dual is the dual of the boundary. 



Definition 3.13. The dual discrete exterior derivative d : Q^{-kK) — )• is defined 

by duality to the boundary operator d : Ck+i{*iK;Z) — )• Ck{K]Z). For a dual discrete form 
S i}^{'kK) and a chain c^+i G Ck+i{MK;7j) we define d by 

{da^,Ck+i) = {a^,dck+i) ■ 

The dual discrete exterior derivative d can be decomposed into the two operators di and 
db, which are respectively dual to di and di, that is 

(da*=, Cfc+i) = (dio^, Cfe+i) + (dbo;^, Cfc+i> = diCk+i) + {a^, ^bCfe+i) • 

The support volumes of a simplex and its dual cell are the same, which suggests that 
there is a natural identification between primal A;-cochains and dual (n — A;)-cochains. 

In the exterior calculus for smooth manifolds, the Hodge star, denoted *, is an isomor- 
phism between the space of A:-forms and (n — A;)-forms. Since the Hodge star operator is 
metric-dependent, in the discrete theory, it is defined as an equality of averages between 
primal and their dual forms [15, 17]. 

Definition 3.14. The discrete Hodge star is a map * : n^^{K) Vq-^{MK) defined by 
its value over simplices and their duals. For a k-simplex , and a discrete k-form a^, 

1,1, !,> I L. 



A^a^^ia") = s-^Aa^a") 



where s is ±1 (see [15]). 



Similarly we can define the discrete Hodge operator on the boundary, that is on an 
(n — l)-dimensional simplicial complex and its dual. It is trivial to show that for a A:-form 
the following holds: * * = (— 

Remark 3.2. The discrete Hodge star can be represented by a matrix (see Section 6). 
According to Definition 3.14, this matrix is diagonal. In the case Whitney forms are used, 
the discrete Hodge operator is sparse but not diagonal in general [17]. 

Next, we define a natural pairing, via so-called primal-dual wedge product, between a 
primal A:-cochain and a dual (n — fe)-cochain. The resulting discrete form is the volume 
form. In order to insure anticommutativity of the primal-dual wedge product, we take the 
following definition. 
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Definition 3.15. Let a'' G ^^i^) be a primal k-form and p"-''' G ^^"''(^i^)- We define 
the discrete primal-dual wedge product A : ^'^{K) x — >■ n^(Vfc(K)) by 

(a*= A^-^V;.) = (l) J^f I (a^a^)(^-^^ia^) 

= (a^(7*^)(^"-^*ia'=) 
= (/^^-'^ A 

where V^k is the n-dimensional support volume obtained by taking the convex hull of the 
simplex and its dual -k^a^ . 

As an illustration, consider the two-dimensional simplicial complex K depicted in 
Figure 1 and its dual *K in Figure 2. Let G ^1{K) and G ^ih^K). The dual cell 

of the primal edge [vk, Vi] is [vm, Vm-i], up to a sign depending on the chosen orientation. 
The primal-dual wedge product of and on the support volume V^^^^^.] = V[vm.Vm-i]^ 
represented by the diamond shaped region generated by [vk,Vi] and [vm,Vm-i], is simply 
a dot product a{[vk,Vi]) ■ P^{[vm,Vm-i])- Now, in order to look at the primal-dual wedge 
product on the boundary, let 7° G and f}^ G fl^{d{-kK)). For instance, fj^ can be a 

restriction of on the bomidary The primal-dual wedge product (7^ A fj^, Vuf.) = 

l^ivk) • fj^{[vp,Vp-i]). The volume for V^^. = V[fj^^^^_-^j is simply the measure of the cell 

[Vp,Vp-l]. 

Here we note the advantage of employing circumcentric with respect to barycentric 
dual since one needs to store only volume information about primal and dual cells, and 

not about the primal-dual convex hulls. 

Definition 3.16. Given two primal discrete k-forms, Q^,j3^ G Q,^{K), their discrete 
inner product, {a^,(3^)d is given by 

= («^a^•)(*/3^*ia^). 

The proposed definition of the dual boundary operator assures the validity of the 
summation by parts relation that parallels the integration by parts formula for smooth 
differential forms. 

Proposition 3.1. Let K be an oriented well-centered simplicial complex. Given a primal 
{k — l)-form a^~^ and a dual (n — k)-discrete form , then 

{da^-^ A K) + (-l)'=-i(a^-^ A d/3'^-^ K) = {a^-^ A dK) , 

where in the boundary pairing a!^~^ is a primal {k — l)-form on dK, while is a dual 

(n — k)-cochain taken on the boundary dual -k{dK). 
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Proof. We have 
and 

(a'^-^ Ad^"-^K) = ^(q!'=-\ (7*^-1) (d^"-^*i(7'=-^) = ^(a'=-^(7*^-^)(/3"-^a(*i(T'=-^)) 

Inducing the orientation of the dual such that s^k = s^k-i = (—1)*^ completes the proof. 

□ 

Remark 3.3. Decomposing the dual form, into the internal and the boundary part 

as B'^~^ = \ J*'"^ I ^'"^ and decomposing the dual exterior derivative in 

1 ySb G fi^-'^Ki^) on d{i.K) 
the same manner, the summation by parts formula can be written as 

(da'^-i AA,i^) + (-l)'^-i(a'^-^A(diA + db/3b),i^) = {a''-' A P^,dK) . (3.1) 

In the standard literature of discrete exterior calculus, the codifferential operator is 
adjoint to the discrete exterior derivative, with respect to the inner products of discrete 
forms [15]. According to Proposition 3.1, this is not the case since on the right a term 
corresponding to the primal-dual pairing on the boundary appears. As the subsequent 
section demonstrates, this term is precisely responsible for the inclusion of the boundary 
variables in the discretized Stokes-Dirac structure. 



4 Dirac Structures on a Simplicial Complex 

The Stokes-Dirac structure, which captures a differential symmetry of the Hamiltonian 
field equations, as presented in [26], is metric-independent. The essence of its construction 
lies in the antisymmetry of the wedge product and the Stokes theorem. In a discrete 
framework, the primal-primal wedge product [15] inherits a number of important proper- 
ties of the cup product [25], such as bilincarity, anti-commutativity and naturality under 
pullback; however, it is degenerate and thus unsuitable for defining a Dirac structure. This 
motivates a formulation of a Dirac structure on a simplicial complex and its dual. We 
introduce Dirac structures with respect to the bilinear pairing between primal and duals 
forms on the underlying discrete manifold. We call these Dirac structures simplicial Dirac 
structures. 
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In the discrete setting, the smooth manifold M is replaced by an n-dimensional well- 
centered oriented manifold-like simplicial complex K. The flow and the effort spaces will 
be the spaces of complementary primal and dual forms. The elements of these two spaces 
are paired via the discrete primal-dual wedge product. Since the Stokes-Dirac structure 
V expresses the coupling between fp and e^, also fg and Cp, via the exterior derivative, 
whose discrete analogue maps primal into primal and dual into dual cochains, the flow 
space cannot be entirely built on a primal simplicial complex and the effort space on a 
dual cell complex, or vice versa. Instead, the flow and the effort spaces will be mixed 
spaces of the primal and dual cochains. One of the two possible choices is 

j-^^g = ^i%mK) X al{K) X aY'idiK)) 

and 

The primal-dual wedge product ensures a bijcctivc relation between the primal and 
dual forms, between the flows and efforts. A natural discrete mirror of the bilinear form 

p,q p,g 



(2.1) is a symmetric pairing on the product space J^ig^ ^va defined by 



((( /p , fl fl 4, e\, el), {pp, fl fl el, e^, 4)))^ 

" ' ' ' . s 

= {e\^^l^e\^^l^el^^l-re\^j\,K)^{ei^fl^el^jl,^K). 

A discrete analogue of the Stokes-Dirac structure is the finite-dimensional Dirac structure 
constructed in the following theorem. 

Theorem 4.1. Given linear spaces Tpg and £pg, and the bilinear form {{,))d- The linear 
subspace C Tp g x £p g defined by 

^d = { (/p) fqi fb: ^p, Cg, 65) G -^p.g X S^ g [ 

A = {-l)Pep\gK} 

is a Dirac structure with respect to the pairing {{,))d ■ 

Proof. In order to show that C V^, let (/^ , /j^, Cp, e^, e^) € V^, and consider any 
ifpJgJb^^'^l^^b) e '^d- Substituting (4.2) into (4.1) yields 

{{-ir-^'el A (die^ + dbeg) + e] A del + {-^^'''4 ^ i^^^l + ^^^D +4^ ^e^, K) 

+{-lf{elAel + elAel,dK). 

By the anticommutativity of the primal-dual wedge product on K 

{elAdel,K) = i-iy(P-'HdelAel,K) 
{elAdel,K) = {-iy^P-'\delAelK), 
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and on the boundary dK 

{eUeldK) = (-l)(f-i)('?-i)(e2 Aei,9i^) 
(elAeldK) = (-l)(f-«('?-i)(ei A e^^ , 

the expression (4.3) can be rewritten as 

(_1).(P-I)(de^ A + {-ir-^el A (d^^ + d^el) , K) 
+(_1)9(P-I)(dei Ael + i-ir-^el A (d^e^ + dbe^) , K) 
+(_1)P+(P-1)(<?-1) (gi A el + el A e^, dK) . 

According to the discrete summation by parts formula (3.1), the following holds 

{del ^ + (- A [die] + d^el) , = (e^ A , OK) 
{del ^ el + {-ir-Pel A {d^el + dbe^) , K) = {el A el dK) . 

Hence, (4.3) is equal to 0, and thus C T)^. 



Since dim J-'p ^ = dimiSp ^ = dimP^, and ((, ))d is a non-degenerate form, Vd = T?^. □ 



Remark 4.1. As with the continuous setting, the simplicial Dirac structure is algebraically 
compositional. Since the simplicial Dirac structure is a finite- dimensional constant 
Dirac structure, it is integrable. 

The other possible discrete analogue of the Stokes-Dirac structure is defined on the 
spaces 

= nr^i^K) X nr'{K) X nr\dK) . 

A natural discrete mirror of the bilinear form (2.1) in this case is a symmetric pairing on 
the product space g x Sp ^ defined by 

1 ^1 ^1 ;1 „1 A\ /J-2 f2 f2 ~2 „2 „2\ 



((( fpi fqi fbi ^g' ^fe)' (/p 1 fqi fbi ^b))) d 



= {elAf^ + elAf^ + elAfl+elAflK) + {elAfi + elAfl,dK). 
Theorem 4.2. The linear space defined by 

'i>d = { ifp, fq, fb, ep, Cy, Cb) G Jp^^y X £^ I 



p,q\ 



;:)Kr"n(::)K:)/" 

efe = (-l)%|aif} 

is a Dirac structure with respect to the bilinear pairing {{, )) j. 

In the following section, the simplicial Dirac structures (4.2) and (4.4) will be used 
as terminus a quo for the geometric formulation of spatially discrete port-Hamiltonian 
systems. 
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5 Port-Hamiltonian Dynamics on a Simplicial Complex 



In the continuous theory, a distributed-parameter port-Hamiltonian system is defined with 
respect to the Stokes-Dirac structure (2.2) by imposing constitutive relations. On the other 
hand, in the discrete framework one can define an open Hamiltonian system with respect to 
the simpHcial Dirac structure or the simphcial structure Dd- The choice of the structure 
lias immediate consequence on the open dynamics since it restricts the choice of freely 
chosen boundary efforts or flows. Firstly, we define dynamics with respect to the structure 
(4.2) and (4.4). Then, in the manner of finite-dimensional port-Hamiltonian systems, we 
include energy dissipation by terminating some of the ports by resistive elements. 

5.1 Port-Hamiltonian systems 

Let a function 7i : fl^{^iK) x Q'^{K) — )• M stand for the Hamiltonian (dtp, aq) i-)- 'H{dip, Og), 
with OLp G Vl^Jy-kiK) and Uq G r2^(K). The value of the Hamiltonian after arbitrary 
variations of dtp and aq for 5ap G ^^J^-kiK) and 5aq G f2^(if), respectively, can, by Taylor 
expansion, be expressed as 

dT-L dT-L 

7i{ap + 5ap, aq + 5aq) = 7i{ap, aq) + {-^ A Sckp + A 6aq, K) 

-|- higher order terms in 5ap,5aq. 

Here, it is important to emphasize that the variations 5ap, Saq are not restricted to vanish 

on the boundary. 

A time derivative of Ti along an arbitrary trajectory t {ap{t) , aq{t)) G 0^(*iK) x 
n'^aiK), t G M, is 

d^,,„ , ,dH da„ 8% duq n\ 
^«K,a,) = {g^A^ + ^A^.i^). (5.2) 

The relation between the simplicial Dirac structure (4.2) and time derivatives of the 
variables are 

t d&p daq r,\ 

^^-^^ 

while the coenergy variables are set 

m . m .... 

ep = eg = — . (5.4) 

dap oaq 

This allows us to define the spatially discrete, and thus finite-dimensional, port- 
Hamiltonian system on a simplicial complex K (and its dual -kK) by 

/ 31-1 \ 

{-ird,\ m; \ , , ,,./db' 

(5.5) 




d M#ii+(-^Hor^' 



d&p 



dK 
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where r = pq + 1. 

It immediately follows that ^ = (^6 A fb, dK), enuneiating a fundamental property of 
the system: the increase in the energy on the domain \K\ is equal to the power supplied 
to the system through the boundary dK and d{'kK). Due to its structural properties, the 
system (5.5) can be called a spatially-discrete time-continuous boundary control system 
with Cfe being the boundary control input and being the output. 

An alternative formulation of a spatially discrete port-Hamiltonian system is given 
in terms of the simplicial Dirac structure (4.4). We start with the Hamiltonian function 
{ap,aq) I—)- Ti{ap,aq), where ap G Q^{K) and G r2^(*iif). In a similar manner as in 
deriving (5.5), we introduce the port-Hamiltonian system 

(t)Kr-r)(|)^U)^'' 

In contrast to (5.5), in the case of the formulation (5.6), the boundary flows can 
be considered to be freely chosen, while the boundary efforts ej, are determined by the 
dynamics. 



5.2 Modelling dissipation 



Incorporation of dissipation parallels the continuous case and for the present moment we 
shall consider only dissipation modeled by port termination. As an illustration, consider 



a mapping : ^liK) — )• '^{*iK) that satisfies 



Furthermore, let = R* with R being a positive real constant. 

In case of the simplicial Dirac structure (4.2), introduce the relation 

= -(-l)«("-?)^d(/,) = -(-l)«("-«)i? * /, , 

as well as associate to every primal {n — ^>)-cell an energy storage effort variable and to 
every dual ^?-cell a sign consistent energy flow leading to 

with H. being a total stored energy. 

This leads to relaxation dynamics of a diffusion process 

^ = (-l)^-^it:di*d^ + (-l)^"'d5e,. 



19 



with 

dH ,dn dx 

= ^ ((-l)^"'^di * + (-irdbe,) , K) 

Let fp = *fp = — fp,x G Q^{K), that is p = n and g = 1. As the stored energy 
take n = l{x A*x,K). Then 

dx 

— = R*di*dx = -RSdx + (-l)"dbe6 = -RAx + (-l)"dbe6 , 
at 

where A is the Laplace operator, A : ^'^{K) -> Q,'^{K). One needs to be careful here with 
the minus sign since by the chosen convention Ax = — divgradx [1]. The boundary flow 
is fb = i-'^)"x\oK- 

Remark 5.1. Consider a diffusion process on a one- dimensional simplicial complex K 

with a dual -kK that is also a one- dimensional simplicial complex (the spatial domain is 
identical to the domain of the telegraph equations, confer to Figure 5). This means that 
n = q = 1 and p = 0. The operator Ra is a positive definite operator that maps the set of 
primal edges into the set of dual nods. The resulting dissipative port-Hamiltonian system 
is 

dx 

— = diRd{dx) 

for Cf, = 0, which conduces to the standard compartmental model. This can he extended to 
structure-preserving discretization of reaction- diffusion systems as hinted in [31]. 



6 Matrix Representations for Linear Port-Hamiltonian Sys- 
tems on a Simplicial Complex and Error Analysis 

Discrete exterior calculus can be implemented using the formalism of linear algebra. All 
discrete fc-forms can be stored into a vector with entries assuming the values that those 
forms take on the ordered set of A;-simplices. The boundary operator is a linear mapping 
from the space of fc-simplices to the space of {k — l)-simplices and can be represented by a 
sparse matrix containing only ±1 elements, while the exterior derivative is its transpose. 
There is a number of different Hodge star implementations, hwl the so-called mass-lumped 
is the simplest, with the Hodge star being a diagonal matrix. 



6.1 Matrix representations of linear operators 

Any discrete differential form G ^^{K) is uniquely characterized by its coefficient vector 
a. G A'^, where A'^ = M^^^, Nj^ = dim^^[K) is the number of /c-simplices. Similarly, for 
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a, Pi E '^{mK) the vector representation is Pi G A„_fe = R^''. Representing discrete 
forms by their coefficient vectors induces a matrix representation for hnear operators (see 
e.g. [9, 17]). 

The exterior derivative d : ^^{K) — )■ 17^^^ (K) is represented by a matrix D'' G 
-^Nk+ixNk ^ which is the transpose of the incidence matrix of /c-faces and {k + l)-faces 
of the primal mesh [9, 10]. The discrete derivative d : Q^{i^K) — > U'^'^^{i^iK) in the 
matrix notation is the transpose of the incidence matrix of the dual mesh denoted by 

f) (z M^fc+ix^fc, which can be, as we shall soon show, decomposed as -D = (Dj : -Db)^ with 
Di and Dh being matrix representations of dj and db, respectively. 

The exterior product A : n'^{K) x n'^-''{iciK) n2{Vk{K)) between a G O^(i^) and 
/5 e can be written as 

(a AP,K)= a'Wj^-''P = (-l)'=("-*^);3'^W?^_fea 
= (-l)*^("-*=)(/3Aa,if), 

where W^'^ ,W^_^ G R^VfexAffe ^n-fc ^ ^_^Y^n-k) (^^k_^y _ 

A crucial ingredient for supplying the result of Theorem 4.2 is a discrete summation 
by parts formula (3.1), which in the context of the simplicial Dirac structures can be 
rewritten as 

(dep A eg, K) + {-lf-^{ep A (d^e^ + dfee^), K) = (cp A e^, dK) , 

where Cp G 9!l'^{K), G Q.'l'^ {mK) , G VL'^'^{-k^K) for q = k and p = n- k + 1. 

Representing the discrete forms by the corresponding coefficient vectors ep G A'^, G 
G An-k+i, /q G A*^, G A^ ^ gives rise to the matrix representation 

of (6.1) 

(^D'^-^epY W^-'^eg + (-l)'=-ie-JW^^"_-*^+i (^""'e, + ^T^e^) 

Here the matrix T*^"^ G M^^'^^fc ^^k-i jg ^ trace operator of (A;— l)-forms on the boundary 
of the primal simplicial complex K. 
After regrouping we have 

Choosing W;^-'' = In,, Wj^Zi^^ = In,., and = /.imA^-i ^P^i^s the well- 

known relation [10] 

D'^-^ = {-if (o^'-^y , 

while interestingly enough the dual boundary operator is a dual of the primal trace oper- 
ator 
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The discrete Hodge operator * : ^'^{K) — >■ ^{-kiK) has the following matrix repre- 
sentation [9, 10] 

M^d = P with M'^ G M^'=^^'= for a e A'' , P e An-k , 
while the discrete Hodge star from A^-k to can be described by 

where M'^~^M^ = Ij^^. As in the continuous theory, the discrete Hodge oper- 

ators are invertible. The norm of a G A'^ induced by the discrete Hodge star is 

\\a\\l^=a'M^a = \\{M'')'2a\\ . 

6.2 Representation of simplicial Dirac structures 

Consider the simplicial Dirac structure (4.2). The effort and flow space are 

J-;f,,A = A,, X A" X A^^ 
and the bilinear form (4.1) on the space J^p^qA x Sp gA is 

IKfl fl fl -1 ^1 ^-IWfS fi 7fi pi ^2n« 

\\\J p 1 J q 1 Jh 1 ) 1 \ J p 1 J q 1 Jh 1 ) II a 

+ {4rK,7-qf! + {4YK7-qfl (6.1) 

+ (-l)("-^)("-^)((e1)V? + (e1)Vl) , 

where we took Wl_p = In,, Wl_^ = (-l)'^("-«)/^^, W^;-^,^ = (-l)("-f)(-<?)7^,_^. 
The matrix representation of the simplicial Dirac structure (4.2) is 



fp ] 




f 


(_l)P9+l£)"-9 ^_ 


fq 


H 







fb J 




I ( — 1)PT"~*' 











(_1)5(P+1)+1 




H 











I (_1)PJ'"-P 






^6 






(6.2) 





/ V efe 
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6.3 Error analysis 

In this section we consider the spatial discretization of a hnear distributed-parameter 
port-Hamiltonian system of the form 



del 

de: 



c 



^ = (6.3) 

on an ra-dimensional polytope \K\. The operators *p and *q are the Hodge stars spawned 
by Riemannian metrics. 

Note that all continuous (spatially undiscretized) quantities are labeled by a superscript 
c, for example, and are the continuous efforts. The approach to convergence analysis 
we take here is that of [18]. 

The discrete analogue of (6.3) defined with respect to the simplicial Dirac structure 
(4.2) is 



-Mqkq = Z)"-f Cp (6-4) 
A = (-l)fT"-%, 

where Mp G M^p^^p and Mq G M^9^^9 are diagonal Hodge matrices. A dot over a variable 
denotes the time derivative. 

Integrate the first equation over dual p-cells and the second over primal g-faces to 
obtain 

-Mptp -kp= (-1)*"?+^ {Dr% + i^r'e":) 

where and are integral forms on the primal mesh and its circumcentric dual, while 
Rp and Rq are time derivatives of the residues of the Hodge operator approximations given 

by 



*qeq = Mq^iCq^i + Rq^i , 



with subscripts k and I acting as selectors for vector components. 

Define discrete energy errors as Sep = — and Seg = — e^, and the output error 

as Sfb = - fb. 
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Subtracting (6.4) from (6.5) leads to 

-Mp5kp - 4 = (-l)P^+i [Dr^Seq + D^'^-Hel) 



-Mg6eq -Rg = D^'-Pdcp 



(6.6) 



since Seb = el-eb = {J^^n-, - {eb,*(^b~^))^^^-i^Q(^^K) = 0. 

Multiplying the first equation in (6.6) by Sep and the second by Scq gives 

~{5ep,Mp6kp) - {5ep,Iip) = {-If'^+^Sep^D'^-Heg) 
-{Seq,MgSkg) - {Seg,Rq) = {Seg,D''-PSep) . 

Then we have 

-{Sep,MpSep) - {Sep,Rp) - {deq,Mqdeq) - {Seq,Rq) 
= {-l)Pi+\dep,D^-''Seq) + {Seq,D^-PSep) =0. 

That is _ . . 

{Sep, MpSep) + {Seq, MqSeq) = -{Sep, Rp) - {SSqRq) . (6.7) 

Integration of (6.7) from to yields 

\\\Mlsep{tf)f + \\\Mlseq{tf)f = - j\sep{r) , kpir)) + (5e,(r), i?,(r))dr 

< /*' ||4(T)||||(5e- (r)|| + \\iiq{T)\\\\5eq{T)\\dT . 
JO 

Let t* be such that 

WMp'Sipinf + WMiSeqinf = max ||Mj|| ||5e- (t)|| + ||m| ||||<5e,|| , 



then 



2 / 



\\M;^5ep{t*)\\ + WMq'SeqinW ) < 2 ( ||M|5e- (f )r + ||M|5e,(r)||' 



< 4 / ' ||4(r)||||,5e- (r)|| + \\iiq{r)\\\\Seq{r)\\dr 

Jo 

< 4^*' (^\\Mp'Sep{t)\\ + \\MiSeq{t)\\^ (^\\MP lip{T)\\ + \\MP iiq{r)\\'j dr . 
It follows that 

\\MlSep{f)\\ + \\Ml5eq{t*)\\ < 4 ||iQ7^i?p(r)|| + \\Mq'^ iiq{T)\\dT . 

Jo 

Thus 

\\Sep{f)\\ + \\Seq{f)\\ < 4 (||il^^||oo + ||M,"^||oo) J^' \\MP lip{T)\\ + ||M,"^4(r)||dr . 
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Estimation of the residues Rp and Rq can be conducted by employing Bramble-Hilbert 
techniques in the case of a weak formulation, or using a Taylor's expansion of the efforts 
under the standard smoothness assumptions [17]. For the results on the estimates of 
the Hodge star in one, two and three dimension the reader is invited to consult [17] and 
references therein. 

7 Physical Examples 

In this section we formulate discrete analogues of distributed-parameter port-Hamiltonian 
systems on a three-, two-, and one-dimensional manifold. 

7.1 Mcixwell's equations 

Let K hea well-centered 3-dimensional manifold-like simplicial complex with circumcentric 
dual -kK, endowed with a discrete Riemannian metric. Mirroring the continuous case [26], 
we formulate the discrete Maxwell's equations in terms of discrete differential forms, and 
then we demonstrate that the underpinning differential/gauge structure is preserved. 

The energy variables are chosen such that they live on the discrete manifolds that are 
dual to one another. For instance, we choose the magnetic (field) induction 2-form to be 
defined on the primal simplicial complex K as aq = B e ^^^(i^) and the electric induction 
2-form ckp = D e Q^{~kiK). This means that B and D do not reside at the same discrete 
locations, but rather at separate faces of staggered lattices. 

Remark 7.1. In the case of a spatio-temporal discretization based on the asynchronous 

variational integrator scheme, as proposed in [33], the electric and magnetic induction are 
also defined at different time locations leading to improved numeric performance (for more 
details refer to [33]). 

The coenergy variables are chosen coherently as implied by the choice of the energy 
variables such that the discrete Maxwell's equations fit the simplicial Dirac structure (4.2) 
for n = 3, p = g = 2. This entails that the magnetic field intensity eq = H e Q,^{-kiK) and 
the electric intensity Cp = E e Cl^{K), as such, are related to the energy variables via 

D = *eE 
B = *iiH , 

where e and ji denote the constant electric and magnetic permittivity, respectively. 
The corresponding simplicial Dirac structure is 




The Hamiltonian isT-L = ^{E A D + H A B, K) , or expressed only in terms of the primal 
forms QS'H = \{E A*eE+^*B AB,K). 
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Under the assumption that there is no current in the medium, the spatially dis- 
cretized Maxwell's equations with respect to the simplicial Dirac structure (7.1) in the 
port-Hamiltonian form are given by 




-di 
d 




m 

dD 

ak 

dB 




(7.2) 



h = 



m 

dt) QK 



The readily proved energy balance is ^ = (cb A /{,, dK). Incorporating a nonzero current 
density into the discrete Maxwell's equations is straightforward as in the continuous case. 

7.2 Two-dimensional wave equation 

In order to demonstrate practically that we do not face a problem of interconnection 
of the elementary Dirac structures encountered in the mixed finite element method, as 
reported by [37] (see pages 183-196), we consider the simplicial Dirac structure behind 
the discretized two-dimensional wave equation. The normalized wave equation is given by 



where ^ is a smooth 0-form on a compact surface M C with a closed boundary, and 
A is the Laplace operator. This equation, together with nonzero energy flow, can be 
formulated as a port-Hamiltonian system with boundary port variables [34, 12]. 

The energy variables of the discretized system are chosen as follows: the kinetic mo- 
mentum is a dual 2-form whose time derivative is set to be /p, the elastic strain is a 
primal 1-form with time derivative corresponding to /q, the cocnergy variables are a pri- 
mal 0-form Cp and a dual 1-form e^. Such a formulation of the discrete wave equation 
is consonant with the simplicial Dirac structure (4.2) for the case when p = n = 2 and 
q = I. We shall, nevertheless, practically confirm the arguments of Theorem 4.2 in a 
simple low-dimensional model. 

Consider a ring of counterclockwise oriented triangles that could be, say, obtained by 
a very coarse discretization of a disk. The dual of the central vertex vq is its Voronoi 
region, while the duals of the boundary vertices are the convex boundary pentagons. The 
orientation of the primal edges is chosen as indicated in Figure 4. The orientation of the 
dual edges is induced such that the basis of the primal and dual cells combined give the 
orientation of the embedding space that, in our case, has been given by the right-hand 
rule (for more on orientation see pages 11-22 of [15]). 

It suffices to check the power conserving property of the founding Dirac structure. We 
need to show that 



{epAfp + egAfg,K) + {ebAfb,dK)=0. 



This is equivalent to the validity of the following relation 



(dep ACq + CpA [diCq + dbC^), K) = (cp A e^, dK) . 
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Figure 4: A simplicial complex K consists of five triangles arranged into a pentagon. The 
dual edges introduced by subdivision are shown dotted. 



We calculate 



{epiv2) - ep{vi)) eq{[vi,VG]) + {ep^vs) - ep{v2)) eq{[v2,V7]) 
+ (ep(v4) - ep{v3)) eq{[v3,vs]) + (ep(v5) - ep(-U4)) eg([t}4, vg]) 
+ (ep(vi) - ep(u5)) eq{[v5,vio]) + (ep(wo) - ep(vi)) eg([i)6, t)io]) 
+ (ep(^^o) - ep{v2)) eq{[v7,V6]) + (ep(vo) - Cpivs)) eq{[vs,V7]) 
+ (ep(vo) - ep(v4)) eq{[vQ,vs\) + (ep(vo) - epiv^)) eq{[vio, vg]) 



and 



(cp A (diCg + dbefe),ir) = ^ (cp, (J°) (die, + dbCfe, 

= E (ep,^°>((e9,5i(*a0)) + (e,,ab(*O» 



= ep(t;i)(^eg([ui,'06]) + eq{[vQ,vio]) - eq{[vr„vio]) + eb{[v5,vi]) 
+ ep{v2)(^eq{[v2,V7]) + eq{[v7,V6]) - e<j( [wi , We] ) + eb{[vi,V2]) 
+ epiv3)(^eq{[v3,vs]) + eq{[vs,V7]) -eg{[v2,V7]) + ebi[v2,V3]) 
+ ep{v4)(^eq{[v4^,V9]) + eq{[vg,vs]) - eq{[v3,Vs]) + efe([z)3, 7)4]) 
+ ep(w5) (eg([%, t'lo]) + eq{[vio, Vg]) - eq{[v4,, Vg]) + ebi[v4, %]) 
+ ep{vo) - eq{[v7, vq]) - eq{[vs, V7]) - eq{[vg, vs]) - eq{[vio, Vg]) 
- eqi[ve,vio])) . 
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After summation of the last two relations, all terms, except those associated with the 
primal and dual boundary, cancel out, leading to 

(dcp A eq,K) + (ep A (dieg + d}yeb),K) = ep(z;i)efe([{)5, ui]) + ep{v2)eb{[vi,V2]) 

+ ep{v?,)h{[v2, vz]) + ep('U4)e6([i)3, 'O4]) (7.3) 
+ ep(t;5)e6(['04,'05]) . 

This confirms that the boundary terms genuinely live on the boundary of \K\. 
7.3 Telegraph equations 

We consider an ideal lossless transmission line on a 1-dimensional simplicial complex. The 
energy variables are the charge density q G Q}^{K), and the flux density (j) G Vl\{-kK), hence 
p = q = 1. The Hamiltonian representing the total energy stored in the transmission line 
with discrete distributed capacitance C and discrete distributed inductance L is 

with co-energy variables: ep = ^ = *^ = F representing voltages and ^5 = ^ = *^ = / 
currents. 

Selecting fp = and fq = leads to the port- Hamiltonian formulation of the 
telegraph equations 




In the case we wanted to have the electrical current as the input, the charge and the flux 
density would be defined on the dual mesh and the primal mesh, respectively. Instead 
of the port-Hamiltonian system in the form (5.5), that is (7.3), the discretized telegraph 
equations would be in the form (5.5). The free boundary variable is always defined on the 
boundary of the dual cell complex. 

Note that the structure (7.3) is in fact a Poisson structure on the state space ^1^{K) x 
0^(*i<r). This will become obvious when we present this structure in a matrix represen- 
tation. Before that, it is illustrative to demonstrate how the pairings between primal and 
dual forms can be rather easily calculated. 

Using the notation from Figure 5, we have 

{deqAep,K) = ^ {deq,a^){ep,*a^) = ^ {eq,da^){ep,ica^) 

= [eq{v2) - eq{vo)]ep{vi) + [eq{v4) - eq{v2)]ep{v3) + ... 

+ [eq{v2n-2) " eg{v2n-4:)]epiv2n-3) 
^p{V2n-l) 

= -eq{vo)ep{vi) - eq{v2)[ep{v3) - ep{vi)] 

- eq{v4)[ep{v5) - ep{v3)] 

- eq{v2n-2)[ep{v2n~l) - ep{v2n-3)] + eq{v2n)ep{v2n-l) 
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^0 V2 f4 



V2n-2 ^2n 



Figure 5: The primal 1-dimensional simplicial complex K with even nodes indices and its 
dual -kK with odd indices, both with conventional orientation of one simplices (from the 
node with a lower- index to the higher-index node). By construction, the nodes % and 
V2n are added to the boundary as previously explained to insure that the boundary of the 
dual is the dual of the boundary, i.e., d{-kK) = •k{dK). 

and 

{eq A (diCp + dy,fb),K) = ^ (e,, o-°)(diep + db/f„*o-°) 

= E (e„cT0)((ep,ai(*O) + (A,c^(*a0))) 

= eqivo)[ep{vi) - fb{vo)] + eq{v2)[ep{v3) - ep{vi)] 
+ eq{v4,)[ep{v5) - ep{v3)] + ... 

+ eq{v2n-2)[epiv2n~l) - ep{v2n-3)] 
+ eg(v2n)[/fe('52n) " ep(l'2„-l)] • 

The arguments of the discrete Stokes theorem (3.1) are trivially verified 
{deq A Cp, K) + {eg A (djCp + dkfh),K) = {cq A A, dK) 

= -eq{vo)fb{vo) + eq{v2n)fh{v2n) 

showing the power-preserving property of the simplicial Dirac structure (7.3) which implies 
that for any (cp, /p, Cg, fq, Cf,, /;,) in the simplicial structure (7.3) the following holds 

(ep A /p, K) + {eq A fq, K) + (e^ A A, OK) = Q . 

The energy balance for the transmission line thus is 
d'H 

= (eb A /b, dK) = eb{v2n) fb{v2n) - eb{vo)fb{vo) , (7.4) 

which demonstrates that the boundary objects genuinely live on the boundary dK. 

Matrix representation. A differential form Cg G Q^{K) is uniquely characterized by its 
coefficient vector Cq G W"-^^ since dimVl^^K) = n + I, similarly Cp, fp G M", fq G M'^"*"^, 
Sb, fb G M?. The exterior derivative d : ^^^(K) — )• Q\{K) is represented by a matrix 
D G M"^("+^), which is the transpose of the incidence matrix of the primal mesh [9, 10]. 
The discrete derivative di : Q^{-kiK) — t- Q\{-kK) in the matrix notation is Di = —D^, 
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and db : ^^^{-^bK) U^{i^K) is represented by Di,, which is the transpose of the trace 
operator. For the simphcial complex in Figure 5, we have 



D 



/ -1 10 
-1 1 

\ GO 

-10 0- 
0- 



\ 




-1 I J 



(7.5) 




1 



Implementing the primal-dual wedge product as a scalar multiplication of the coeffi- 
cient vectors and by taking for convenience e^, = (efe(fo), — efe(f2n))^; the simplicial Dirac 
structure of (7.3) can be represented by 



(7.6) 






-Dl 



Expressing flows in terms of efforts and choosing L 
matrix formulation of (7.3) is 

-MpCp = Deq 

-Mqkq = -Djep + Dbh 



C 



1 for convenience, the 



(7.7) 



where Cp e Aq, Sg G A^, Mp = diag(/ii, hs, . . . , /i2n-i) € M"^" and Mg = diag^^o, ^2, • • ■ , ^2n) € 

^(n+l)x(n+l)^ with hi = \[vo,V2]\, H = \[v2,Va]\, /i2n-l = | ^2^-2, ^^2n] | and ho = 

\[vQ,Vi]\, h2 = \[V1,V3\1 ■■■,h2n = \[v2n-l, V2n\\- 

Error analysis. The time derivatives of the Rp components are 



Rp,l — h2i-iep{v2i-i) 



'[V21-1,V21 

The Taylor's expansion of around V21-1 is 



*e'p = h2i-iep{v2i-i) 



V21 



Q-c Ft ( " ^2 

eliz) = e;{v2i-i) + -Q^iv2i-i)iz - V21-1) + ^(^2^-1) ^" ^"^^ + 0{z') 



Thus 



h2i 



/V21 I c^2^c 

e;iz)dz = el{v2i-i)h2i-i + 3^(*2;-i) ( ^ ) + Oih\ 



2i) ) 



and 



1. 1 d'^e'' 



21 



+ o(4). 
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Likewise, 



Rq,k+1 = h2keqiv2k) - = /i2iteq(u2jk) - / eqiz)dz , 

where k = 0,1, . . . ,n and V-i = vq. In the similar fashion, we obtain 



Rq,k+i = ^^(^2fc) ( ^ j +0{hlk) for k = l,...,n-l 



^ 1 ^e"" 

R,.i = i^^{vo)hl + O{hl) 
% 1 de^ 

Rq,n+1 = -2 (^2n)/iL + 0(^L) • 

It follows that there exist Ki , K2 G M such that 



f f 

\\Sep{f)\\ + \\deg{t*)\\<hK2 / llr^ 

Jo 




Hr)\\ + \\^ir)\\]dT + 0{h^) 



(7.8) 



where T is a trace operator, T = D^, and h = mm{hi,hs, . . . , h2n-i, ho, h2, ■ ■ ■ , h2n}- 

Numerical example. To evince how exactly the discrete model relates to the continuous 
one, we take the example from Section 5 in [12]. 

The spatial domain of the transmission system is the line segment M = [0,e — 1]. 
The distributed capacitance and the distributive inductance are z i-> ^^{z) = and 
z L'^iz) = j^, z G M. On the left-hand side a causal input voltage t u{t) 
is assigned, and at the other end the transmission line is terminated by a load of unit 
resistance, meaning *q{t,c — 1) = e — 1). Initial conditions are assumed to be zero, 
i.e. q'^{0, z) = and ^'^(0, z) = for z e Z. The exact solution for the voltage distribution 
is {t, z) ^ e^p{t, z) = u{t - ln(z + 1)), for t > 0. 

Using equidistant division of M and diagonal Hodge operators, the results of numerical 
simulation when the input 6^(0, i) = u{t) = sint, t > 0, are given in Figure 6. The time 
integration technique is Runge-Kutta 4 and the integration step is 0.01. 

All numerical experiments indicate that the discrepancy between the exact value of 
the voltage and the value obtained by numerical simulation is the greatest at the spatial 
point z = e — 1. Thence, in the left-hand side of Figure 6 we show this error as a function 
of time. In all computational experiments, this error, similar to the results in [12], exhibits 
an oscillatory behavior with the amplitude not exceeding the maximum displayed in the 
first period. 

Repeating simulation experiments for uniform grids of different densities indicates that 
the accuracy of the proposed method is 1/n, what comes as no surprise since we worked 
with diagonal Hodge operators, which are of first-order accuracy as shown in (7.8) for the 
system (7.7). 
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Figure 6: On the left, the voltage distribution Cp for n = 10, and on the right, the voltage 
error at the point z = e — 1, that is eb{v2o){t) — sin(t — 1) for t > 0. 

8 Concluding Remarks 

In the framework of discrete exterior calculus, we have established the theoretical founda- 
tion for formulation of time-continuous spatially-discrete port-Hamiltonian systems. The 
staple fiber of our approach is the formulation of the simplicial Dirac structures as dis- 
crete analogues of the Stokes-Dirac structure. These discrete finite-dimensional Dirac 
structures are the foundation for the definition of open finite-dimensional systems with 
Hamiltonian dynamics. Such an approach to discretization transfers the essential topo- 
logical, geometrical, and physical properties from distributed-parameter systems to their 
finite-dimensional analogues. By preserving the Hamiltonian structure, this methodology 
utilizes the analysis and control synthesis for the discretized systems. 

A number of interesting topics and open questions still need to be addressed. Here we 
provide a few miscellaneous reflections and some comments on future work. 

Numerical aspects. The discrete exterior calculus employed in this paper is founded 
on the idea of a simplicial complex and its circumcentric dual. While for some problems 
Delaunay triangulation is desirable since it reduces the maximum aspect of the mesh, 
for others the construction of circumcentric duals might be too expensive (see [29] and 
references therein) . This motivates the development of a discrete calculus on non-simplicial 
complex meshes, such as a general CW complex [16] or a rectangular scheme. Although 
the latter might be inappropriate for geometrically complex objects, a potential advantage 
would be its conceptual simplicity since the circumcentric dual is again a rectangular mesh. 

A major challenge from the numerical analysis standpoint is to offer a careful study 
of the convergence properties of discrete exterior calculus. Furthermore, it would be 
desirable to have higher-order discrete analogues of the smooth geometric operators. This 
primarily pertains to deriving higher-accuracy Hodge star operators, which would possibly 
in return make structure-preserving discretization more competitive even in the domains 
where structure is put aside. A recent article [2] reports some significant initial results 
regarding stability of finite element exterior calculus. The abstract theory is applied to 
linear elliptic partial differential equations with intention to capture the key structure of 
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de Rham cohomology and as such mainly pertains to the vanishing boundary constraints. 
Another related publication [19] extends the framework of [2] to approximate domains. In 
the future, in the context of [2, 19], it would be interesting to study structure-preserving 
discretization of port-Hamiltonian systems in the framework of Hilbert complexes. 

Open discretized systems. The Stokes-Dirac structure has proven to be successful in 
capturing the essential geometry behind many open systems with Hamiltonian dynamics. 
The concept of the Stokes-Dirac structure as presented in the introduction in order to 
accommodate some port-Hamiltonian systems, such as the ideal isentropic fluid, needs 
to be augmented [26]. The main idea behind these modifications remains to be based 
on the Stokes theorem. Prom a structure-preserving discretization point of view, there 
appears not to be any impediments; nonetheless, in order to discuss these questions in a 
systematic manner, a unified theory of open infinite-dimensional Hamiltonian systems is 
needed. The main novelty in discretizing some so formulated general underlying structures 
might concern their integrability. The simplicial Dirac structures formulated in this paper 
are constant Dirac structures and as such they satisfy the usual integrability conditions 
[6, 7, 11]. 

An important application of structure-preserving discretization of port-Hamiltonian 
systems might be in (optimal) control theory, what also prompts a need for time dis- 
cretization. For closed Hamiltonian systems, it is well-known that asynchronous vari- 
ational integrators in general cannot preserve the Hamiltonian exactly; however, these 
integrators, for small time steps, can preserve a nearby Hamiltonian up to exponentially 
small errors [20, 21, 22, 23, 24]. An important issue in this context is to study the effects 
these integrators have on passivity (and losslessness) of open dynamical systems. 

Covariant formulation. It is known that, for instance. Maxwell's equations are also 
consonant with multisymplectic structure since they can be derived from the Hamiltonian 
variational principle [1, 22]. The multysymplectic structure behind Maxwell's equations, 
unlike the Stokes-Dirac structure, is defined, not on a spatial manifold M, but on a 
spacetime manifold X. Here we need to notice that one could define a Stokes-Dirac type 
structure on a pseudo-Riemannian, say Lorentzian, manifold. In Lorentzian spacetime, the 
forms E and B can be combined into a single object, the Faraday 2-form F = E Adt + B. 
The form F can also be expressed in terms of the electromagnetic potential 1-form A 
as F = dA. The Hodge star of F is a dual 2-form G = *F = H A dt — D, known as 
the Maxwell's 2-form. The charge density p and current density J can be combined into 
the source 3-form j = J A dt — p. A well-known relativistically covariant formulation of 
Maxwell's equations [1] is: dF = and dG = j. 

In order to relate this formulation to the port-Hamiltonian framework, define the 
following Stokes-Dirac structure on a Lorentzian manifold X by 

= {ifp, fq, fb, ep, eg, Cb) e 

n'^ix) X n^{x) X n^{dx) x n^{x) x n^{x) x n^{dx)\ (s.i) 

fp = -dCq , fq = dCp , Cb = -eq\dx , fb = ep\dx} ■ 

Similar to the proof of Theorem 1, it is easy to verify that Dc = T^c^-, with respect to 
a natural bilinear form. The dynamics of Maxwell's equations can now be imposed by 
setting fp = F, fq = j, Cp = A, and Cq = G. Furthermore, since d^ = 0, it follows that 
dfp = ddeq = 0. 



33 



A natural choice for discretization of the structure T>c, in the context of discrete 
exterior calculus, would be on a simplicial 4-complex. This would insure a completely 
covariant formulation of discrete Maxwell's equations, similar to Regge's formalism for 
producing simplicial approximations of spacetime in numerical general relativity. The 
relativistic effects in most engineering applications are however negligible, hence for these 
purposes, by choosing a time coordinate, we can split the Lorentzian manifold into 3 + 1 
space, whose discrete analogue is a prismal cell complex. Similar to discretization of 
multisymplectic structures, this would lead to a certain type of asynchronous variational 
integrator. 

An important and challenging avenue for future work is to make an explicit relation 
between multisymplectic and Stokes-Dirac structures, and then to compare their discrete 
analogues. 
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